//===========================================================================
/*!
 *  \file Individual.h
 *
 *  \brief  Parent class of Molecular or Atomic clusters
 */
//===========================================================================
#ifndef _Individual_H
#define _Individual_H

#include "comm/def.h"
#include "comm/matrix_macro.h"

#include <math.h>
#include <iomanip>
using namespace std;
#include <iostream>
#include <vector>
#include <boost/format.hpp>
#include <sstream>				 

using namespace boost;

enum FILE_FORMAT{ F_XYZ, F_GAUOUT, F_GAUOUT_C,F_GAUIN, F_GRO ,F_PDB, F_VASPOUT, F_CAR, F_LATTICE, F_INTERNAL };//!<\enum FILE_FORMAT

enum SPECIES{ ATOMIC=0, WATERCLUSTER=1, LATTICE=2};
	
enum MODEL{ NONE=0, TIP3P=1, TIP4P=2,TIP5P=3 };
enum CORTYPE{ DIRECT = 0, FRAC=1};

enum ATOMICNAME{ UNDEFINED=0, HYDROGEN=1, HELIUM=2, LITHIUM=3, BERYLLIUM=4, BORON=5,
	CARBON=6, NITROGEN=7, OXYGEN=8, FLUORINE=9};
	
typedef struct MAInput1 {
	vector<string> info;
	
	MAInput1(vector<string> info_){ info=info_;}
	
	int getDouble(string key,double &a){
		for(int i=info.size()-2;i>=0;i--)
			if( info[i] == key ){
				a= atof(info[i+1].c_str());
				return 1;
			}
		return 0;
	}
	
	int getInt(string key,int &a){
		for(int i=info.size()-2;i>=0;i--)
			if( info[i] == key ){
				a= atoi(info[i+1].c_str());
				return 1;
			}
		return 0;
	}
	
	int getString(string key,string &a){
		for(int i=info.size()-2;i>=0;i--)
			if( info[i] == key ){
				a= info[i+1];
				return 1;
			}
			return 0;
	}
	
} MAInput1;

typedef struct {
	int t1,t2;
	double dmin,dmax;
} Bond2Type_str;
	
typedef struct Genes{
	public:
		double a[DGR],alpha[DGR]; //lattice parameters
		double x[SIZE_MAX*DGR]; //coordinates
		double force[SIZE_MAX*DGR];
		double freq[SIZE_MAX*DGR];//!< frequecies (cm-1)
		double energy,rmsGrad;
		double USRsig[12];	//for USR
		int iCorType; //!< direct or fractional coordinates
		int Nz[SIZE_MAX];
		int tag; //!< handle of Individual 
		int nAtom;
		int nType[TYPE_MAX];
		int nFreq; //!< number of eigen values or frequencies
} Genes_str;


#define FOR_ATOM(i,t) for(int (i)=0;(i)<nAtom;(i)++) if(Nz[(i)]==(t))
/*!
 *  \class Individual
 *  \brief Parent class of Molecular or Atomic structures
 */
class Individual: public Genes_str{
	public:
		
		Individual();
		virtual void clear();
		
		virtual ~Individual(){};

		//cluster operators
		void RotateX(double theta);//!< rotate around X axis
		void RotateY(double theta);//!< rotate around Y axis
		void RotateZ(double theta);//!< rotate around Z axis
		void RotateAxis(double theta,double axis[]); //!< rotate around an arbitrary axis indicated by a unit vector
		void Translate(double *vec,int direction=1); //!< translate by a constant vector
		void Transform(double D[][DGR]);
		void Transform(double **D);

		//atomic operators
		void RotateX_atom(int i,double theta);
		void RotateY_atom(int i,double theta);
		void RotateZ_atom(int i,double theta);
		void RotateAxis_atom(int i,double theta,double axis[]); //!< rotate a specified atom around an arbitrary axis indicated by a unit vector
		void RotateAxisWithPoint_atom(int i,double theta,double axis[],double point[]);
		
		void Translate_atom(int i,double *vec,int direction=1);		
		void Transform_atom(int i, double D[][DGR]);//!< perform transformation defined by D
		void Transform_atom(int i, double **D);//!< perform transformation defined by D


		//molecular operators
		void Translate_mol(vector<int> &m,double *vec,int direction=1);
		void RotateAxis_mol(vector<int> &m,double theta,double axis[]); //!< rotate a specified atom around an arbitrary axis indicated by a unit vector
		void RotateAxisWithPoint_mol(vector<int> &m,double theta,double axis[],double point[]);
		
		void FixOrigin(); //!< translate to the mass center
		virtual void StdOrient(); //!< rotate to the standard orientation
		virtual void AddStructure(Genes_str &p);
		virtual void Migrate(Genes_str &p);
		void Migrate_atom(Genes_str &p,int from,int to){
			p.Nz[to]=Nz[from];
			p.x[to*DGR+0]=x[from*DGR+0];
			p.x[to*DGR+1]=x[from*DGR+1];
			p.x[to*DGR+2]=x[from*DGR+2];
		}
				
		virtual void Get(Genes_str &p);

		inline double mass(int i){ return Mass[Nz[i]];} //!< return mass of an atom

		void GetMolWithH(int iO,vector<int> &mol){		
			mol.clear(); mol.push_back(iO);
			for(int i=0;i<2;i++)	
				if(Nz[iO]==HBondList[i].t1){ 	
					for(int j=0;j<nAtom;j++) if(Nz[j]==HYDROGEN)
						if(bBond2(iO,j,HBondList[i])) mol.push_back(j);
					break;
				}
		}
		
		//retrieve data operator
		Individual& operator=(const Individual& p);
		double* operator[](int i){ return x+i*DGR;}
		bool operator>(const Individual& right){	return (this->energy > right.energy);	}
		bool operator<(const Individual& right){	return this->energy < right.energy;	}
		bool operator==(const Individual& right){	return this->energy == right.energy;	}
		
		
		
		double CalcUSRsimilar(Genes_str &p);
		bool CalcUSRsig(); //!< calculate USR signature
		double CalcRms(); //!< calcuate RMS based on force[]
		void CalcMassCenter(double rc[]); //!< calculate Mass center of structures
		void CalcInertiaTensor(double I[][DGR]); 
		void CalcPrincipalMoments(double Ip[],double X[][DGR]);//eigen values and eigen vectors of I
		
		int Read(istream &is,int format);
		void Write(ostream &os,int format);
		
		
		virtual void CorrectOrder();//!< correct order of atoms
		virtual bool bCorrect(string &error);//!< check the validality of structures
		//bool bCorrect(){ string error=""; return bCorrect(error);}
		
			
		static const char *Elements[];//!< Name of elements 
		static const char *AtomicNo[];
		static const double Mass[];
		static Bond2Type_str HBondList[];
		static Bond2Type_str bondList[];

		static Bond2Type_str bondOH,bondHO,bondOO,bondFH,bondFF,bondCO,bondCH;
		
		void TIP_convert(int ow,int hw1,int hw2,Individual &water,int imodel);

		//distance and vector operators
		virtual void dR(int i,int j,double *dr){
			for(int l=0;l<DGR;l++)	dr[l]=x[j*DGR+l]-x[i*DGR+l];
		};
		
		void R(int i,double *dr){ for(int l=0;l<DGR;l++) dr[l]=x[i*DGR+l]; }

		void dnR(int i,int j,double *dr){ //!< normalised dR(i,j)
			dR(i,j,dr);	double d=DOTPRODUCT(dr,dr);
			for(int l=0;l<DGR;l++)	dr[l]/=d;
		};
		
		double distance(int i,int j){//!< distance between i and j atoms
			double vec[DGR];	dR(i,j,vec);
			return  sqrt(DOTPRODUCT(vec,vec));
		}
		
		double distance(int i,double *d=NULL){//!< distance between i and a given point
			double vec[DGR]={x[i*DGR],x[i*DGR+1],x[i*DGR+2]};
			if(d!=NULL){ vec[0]-=d[0];vec[1]-=d[1];vec[2]-=d[2];}
			return  sqrt(DOTPRODUCT(vec,vec));
		}
		
		double angle(int ow,int hw1,int hw2){ //!< angle between (hw1,ow,hw2)
			double l1[DGR],l2[DGR]; dR(hw1,ow,l1); dR(hw2,ow,l2);
			return acos(DOTPRODUCT(l1,l2)/sqrt(DOTPRODUCT(l1,l1)*DOTPRODUCT(l2,l2)));
		}
		
		double angle_dgr(int ow,int hw1,int hw2){ //!< angle between (hw1,ow,hw2) in degree
			return angle(ow,hw1,hw2)*180.0/M_PI;
		}
		
		double torsion(int i1,int i2,int i3,int i4); //!< torsion angle of four atoms i1-i2-i3-i4
		
		double torsion_dgr(int i1,int i2,int i3,int i4){ //!< torsion angle of four atoms i1-i2-i3-i4 in degree
			return torsion(i1,i2,i3,i4)*180.0/M_PI;
		}

		static string trim(string s);
		
		bool bBond2(int i1,int i2,Bond2Type_str  &bondType);//!< return true if i1 and i2 are bonded according to the rule defined by bondType
		
		void Remove_atom(int iAtom){
			if(iAtom<nAtom){
				for(int i=iAtom;i<nAtom;i++)	swapAtom(i,i+1);
				Nz[nAtom-1]=0;	nType[Nz[iAtom]]--;		nAtom--;
			}
		}

		void swapAtom(int i,int j){ //!< swapping i and j atom including force
			int tempI,l;
			double tempD;
			tempI=Nz[i];	Nz[i]=Nz[j];	Nz[j]=tempI;
				
			for(l=0;l<DGR;l++){
				tempD=x[i*DGR+l];
				x[i*DGR+l]=x[j*DGR+l];
				x[j*DGR+l]=tempD;

				tempD=force[i*DGR+l];
				force[i*DGR+l]=force[j*DGR+l];
				force[j*DGR+l]=tempD;
			}			
		}
		
		static double Charge(int iType){
			switch(iType){
				case HYDROGEN: return 1;
				case OXYGEN: return -2;
				case FLUORINE: return -1;
				default: return 0;
			}
		}
		
		double TotalCharge(){
			double q=0;
			for(int i=0;i<nAtom;i++) q+=Charge(Nz[i]);
			return q;
		}
		
	protected:
		int getAtomicType(string &s);
		
		//vector< vector<int> > molList;
		//virtual int UpdateMolList(){};
		
		
		virtual int ReadXYZ(istream &is);
		virtual int ReadGro(istream &is);
		virtual int ReadPDB(istream &is);
		virtual int ReadCAR(istream &is){};
		virtual int ReadGauOut(istream &is);
		virtual int ReadGauOut_C(istream &is);
		virtual int ReadVASPOut(istream &is);
		
		virtual void WriteXYZ(ostream &os);		
		virtual void WriteCAR(ostream &os){};
		virtual void WriteGro(ostream &os);
		virtual void WritePDB(ostream &os);
		virtual void WriteGauOut(ostream &os){};		
		virtual void WriteXYZLattice(ostream &os);
		virtual void WriteInternal(ostream &os);
		
		bool bConnectivity(int d[][SIZE_MAX],int n);

		int jacobi33(double a[][3],double d[], double v[][3]);	
};


#endif